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Abstract. We consider the initial value problem for the spheiically symmetric, focusing 
cubic wave equation in three spatial dimensions. We give numerical and analytical evidence for 
the existence of a universal attractor which encompasses both global and blowup solutions. As 
a byproduct we get an explicit description of the critical behavior at the threshold of blowup. 



1. Introduction 

We consider a semilinear wave equation in three spatial dimensions with the focusing cubic 
nonlinearity 

ducj} - A</. - 03 = . (1) 

Heuristically, the dynamics of solutions of this equation can be viewed as a competition 
between the Laplacian which tends to disperse the waves and the nonlinearity which tends to 
concentrate the waves. For small initial data the dispersion "wins" leading to global solutions 
which decay to zero as t ^ oo. For large initial data the dispersive spreading is too weak 
to counterbalance the focusing nonlinearity and solutions blow up in finite time. For each of 
these generic evolutions the leading asymptotic behavior is known: small solutions decay as 
1 /t"^ at timelike infinity while large solutions diverge as \/2/ {T — t) for t approaching 

a blowup time T |I3] I?]. The dichotomy of dispersion and blowup brings up the question of 
what determines the boundary between these two behaviors and, in particular, what is the 
evolution of critical initial data which lie on the boundary. During numerical investigations of 
this question we observed that for a large set of initial data the solutions rapidly converge to a 
universal attractor which is given by a two-parameter family of explicit solutions of equation 
([T]). The aim of this paper is to give analytic and numerical evidence for this behavior (which is 
rather surprising for a conservative wave equation). The description of the critical dynamics, 
which motivated our investigations, emerges as a special case. 

The paper is organized as follows. In section 2 we recall some basic properties of 
solutions of equation ([U and we formulate three conjectures about the existence of the 
attractor Analytic evidence for these conjectures is given in section 3. After a short 
description of our numerical methods and tests in section 4, we present numerical evidence 
for the conjectures in section 5. Finally, in section 6 we make some general remarks. 



Universality of global dynamics for the cubic wave equation 



2 



2. Preliminaries and conjectures 

In this paper we restrict our attention to spherically symmetric solutions </) = (/)(t, r), so 
equation ([U reduces to 

2 

dtt4> — drr<l> drCj) — (j) — . (2) 

r 

This equation is invariant under the following transformations: 

• translation in time Ta by a constant a 

Ta : 0(t, r) ^ (f,{t + a,r), (3) 

• scaling S\ by a positive constant A 

S,:Ht,r)^'-^ (4) 

• conformal inversion / (which is an involution) 

I--<^it^r)^^^^-^,^y (5) 

• reflection (j) —(j). 

Neglecting the Laplacian in ([U and solving the ordinary differential equation dtt(l> — 0^ = 0, 
one gets the one-parameter family of spatially homogeneous solutions 

0b = 7^ . (6) 
This family is a special case of the two-parameter family of solutions ||5| 

4>{a,b){t,r) = — ^ , (7) 

^ ' t + a + b[{t + ay — r^) 

which can be obtained from (|6]l by the action of conformal inversion / followed by the time 
translation Ta- The scaling transformation 5a acting on (t>{a,b) only rescales the parameters 
of the solution without changing its form. Note that the solution (j)[a.b) {ti t) is singular on the 
two-sheeted hyperboloid t = -a - 1/(26) ± ^1/(462) + r^. 

The main result of this paper is the observation that the family (|7]i is an attractor for a 
large set of initial data. More precisely, we have the following conjectures about the forward- 
in-time behavior of solutions of equation (|2j starting from smooth, compactly supported (or 
suitably localized) initial data (by time reflection symmetry analogous conjectures can be 
formulated for the backward-in-time behavior): 

Conjecture 1. For any generic globally regular solution (f){t, r) there exist parameters 
(a, 6) £ R X M+ and k = ±1 such that {(f>{t, r) — k <p(a,b) {t, t)) is bounded for all finite r 
and t — > cx). 

Conjecture 2. For any solution (j){t, r) which blows up at the origin in finite time T, there 
exist parameters (a, 6) G R x R~ and k = ±1 such that (T — t)^'^ {(f>{t, r) ~ k (f>(a,b) {t, 
is bounded for t ^ inside the past light cone of the blowup point. 

Conjecture 3. The borderline between dispersive and blowup solutions, described respectively 
in Conjectures 1 and 2, consists of codimension-one globally regular solutions (j){t, r) for 
which there exist a parameter a e R and k = ±1 such that {(f){t,r) — K'4'(afi){iif)') is 
bounded for all finite r and t ^ co. 
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The parameters (a, h) for which the above assertions hold will be referred to as optimal. 

Remark 1: Conjectures 1 and 2 are refinements of the well-known asymptotic behavior of 
solutions of equation (|2|l: t~'^ decay near timelike infinity in the case of global regularity 
(for small initial data) and \/2/{T — t) growth in the case of blowup (for large initial data), 
respectively. 

Remark 2: Although Conjecture 3 is a special limiting case of Conjecture 1, we state it as a 
separate conjecture to emphasize the critical character of the attractor solution 0(a.o)- Note 
that the slowly decaying global solutions described in Conjecture 3 are not asymptotically 
free, i.e., they do not scatter 

Remark 3: The genericity condition in Conjecture 1 is essential because, as we shall see 
below, there do exist non-generic very rapidly decaying globally regular solutions which do 
not converge to the attractor 4>(^a.b)- 

Remark 4: For very large amplitudes the blowup first occurs on a sphere H] and then 
Conjecture 2 does not hold. 

In the remainder of the paper we give evidence for the above conjectures. The evidence 
is based mainly on numerical simulations, however, before presenting numerics, we will give 
two analytic arguments: one based on linearized stability analysis and one based on an explicit 
solution. 



3. Analytic evidence 

3.1. Linearized stability 

In this section we discuss the linearized stability of the attractor solution (j>(^a,b){ti f)- Since 
linearization commutes with symmetries, we may set a = 6 = without loss of generality. 
We denote the resulting solution by 00, thus 00 = In order to determine the spectrum 

of small perturbations around this forward self-similar solution we will use the symmetry 
under conformal inversions. To this end, consider the region /+(0), the interior of the future 
light cone of the origin (0, 0) of the Minkowski spacetime, and its foliation by spaceUke 
hyperboloids (in this paragraph we follow Christodoulou (|6l) 

t = c + + r2 , (8) 

where c is a positive constant. The hyperboloid ([8]) is asymptotic to 9J+(c), the future light 
cone of the point (c, 0). Let us make the conformal inversion 

/:(t,r)^(i>-)=(^-^,^^) . (9) 

This transformation maps I^{0) (in the original coordinate system) to I^{0) (the interior 
of the past light cone of the origin in the barred coordinate system). In particular, the 
hyperboloids ^ are mapped to spacelike hyperplanes 

i = —c, where c = -7- . (10) 

2c 

Now, consider a global-in-time solution (/)(t, r) inside I^{0). In the barred coordinates 
we have from (|5]l 

0(t,f) - (t2_r2)0(i,r). (11) 

It follows from the above paragraph that the study of the asymptotics of 0(t, r) for t ^ 00 
is equivalent to the study of the asymptotics of (f>{t,f) for f ^ 0^. In particular, there 
is equivalence between linearized perturbations of 00 = \/2/t for t ^ 00 and linearized 
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perturbations of 0o = V2/ {—'t) for < — > 0^. But the latter have already been determined in 
IHIItI. Namely, it has been shown there that the spectrum of smooth linearized perturbations 
about the backward self-similar solution v^/ (— consists of a discrete set of eigenmodes of 
the form (-t)^"Cn(2/)' where y = f/{-t) and Ao = -2,Ai = 0, A„ = n (n > 2). More 
precisely, we have the following eigenmode expansion around (f>o 

54>{t, f) = coi-i)-^ + ci(l - y') + C2i-tfil - + ^^4) ^ o{{-t)') . (12) 

By (fTTl i. this translates into the following eigenmode expansion around (where now 
y = r/t) 

S4>{t,r) = co(l - y^) + cil + c,l ^'y_/y/^ + 0{l/t') . (13) 

The first two eigenmodes in the expansions (fT2] i and (T3[ correspond to the perturbations 
along the symmetry orbits 0(a,&) {t, r) and (f>[a,b) {ti respectively. For example, we have 

d 9 1 

^<A(a,6)(i,^)l(a=0,b=0) ~ 1 - , ^0(a,6) (i, J") I (a=0,b=0) ~ ■ (14) 

The choice of the optimal parameters (a, 6) for the attractor amounts to tuning away the 
coefficients of the symmetry modes cq and ci, hence the rate of convergence to the attractor 
is expected to be governed by the third eigenmode in the expansions (fT2] | and ( fTSI l: 

Ht, r) - 4>ia,b) it, r) ^ for t ^ 0^ , (15) 

and 

r) - </.(,,b) (t, r) ^l/t" for t ^ oo , (16) 

which is consistent with our conjectures. Below we will verify this expectation numerically, 
but first we want to give a simple example which corroborates ( fT6t . 



3.2. The conformal solution 

Equation (|2]i has the following globally regular explicit solution 

2 

(t>conf{t,r) = , (17) 

v/(l + (i-r)2)(l + (t + r)2) 

which we will refer to as the conformal solution. Note that this solution corresponds to time 
symmetric initial data 

0(O,r)--A^, 9t0(O,r)=O. (18) 
1 + 

The conformal solution can be easily found in the framework of conformal compactification 
which maps Minkowski space with the flat metric 77 into the Einstein universe with the 
conformal metric g = il^/y, where the conformal factor 51 = 2[(H-(< — r)^)(l + (t+r)^)]~^/^ 
(S). Under this mapping the cubic wave equation □(/)+(/)^ = transforms to Dg^ — $ + $^ = 
0, where $ = Vl^^cf). The trivial constant solution $ = 1 gives 0co)i/ = ^■ 

By elementary calculation we find that for the conformal solution the optimal parameters 
of the attractor (a, 5) = (— 1/\/2, l/\/2)- More precisely, we have for t ^ 00 

cl>.onf{t, r) - 0(-^.-^)(t, r) = -1 ^i-t^ + 0{t-^) , (19) 
in agreement with Conjecture 1. 
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Figure 1. The future domain of a t = surface is partially depicted on the left panel in 
standard coordinates and entirely represented on the right panel in a Penrose diagram 
Dashed lines are level sets of t, solid lines are hyperboloids shifted in time as given by )20K the 
thick straight lines depict outgoing characteristics indicating the location of the approximate 
wavefront. 



4. The numerical method and tests 



4.1. The hyperboloidal initial value problem for the cubic wave equation 

In order to study the asymptotic behavior of solutions and verify the conjectured convergence 
to the attractor, it is convenient to foliate Minkowski spacetime by hyperboloids. A time 
coordinate r adapted to such a foliation can be written as 



t-^^^+r^. (20) 

The level surfaces of r are standard hyperboloids shifted in the time direction. They have 
constant mean curvature K which is a free parameter. In the following we set K = 3 for 
simplicity. A foliation of Minkowski spacetime by level sets of r is depicted in Figure[T] 

In order to be able to analyze the propagation of the outgoing waves to infinity we 
introduce a compactifying radial coordinate p along the surfaces of our foliation. For the 
regularity of our equations in this compactified setting, we perform a suitable conformal 
rescaling of the metric. The rescaling factor denoted by fl must be a function of to ensure 
regularity at the origin in the conformal manifold. We choose = (1 — /9^)/2 following 
lfm[T2l [T3l. The compactifying coordinate is then chosen according to r = p/fl. 

We express the standard Minkowski metric r/ in the new coordinates (r, p) and rescale it 
with the conformal factor 17^ to obtain 

g = n^rj = -Vi^ dr^ - 2pdTdp + dp^ + p^ da^ , (21) 

where da^ is the standard metric on the unit two-sphere. In these coordinates the zero set of 
the conformal factor corresponds to future null infinity denoted by We rewrite the cubic 
wave equation on Minkowski spacetime in conformally covariant form 

□g$ - ^R[g] $ + ^'"^ = 0, where * = ^ , .9 = ^^V ■ (22) 
The Ricci scalar of the metric g from (l2Tl i appearing in the above equation is given by 
12(l-p^)(3 + p^) 

""^'^^ — {TT^r — • 

With the auxiliary variables, 

2 

ij:=dp<i> and tt := -{dr^ + pdp^), 

1 + P'' 
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Figure 2. Convergence factors in time measured in the L2-riorm. The convergence factor Q 
114,'°™-$" 



is defined by Q ■ 



On the left panel we see that after a short transient 



phase the code converges with 6th order. The dashed curve corresponds to a simulation with 
a 10 times smaller Courant factor and shows 6th order convergence from the beginning. The 
right panel shows results of a long time convergence test and indicates loss of convergence 
after about r = 1000. 



we can rewrite the system (l22l i in first order symmetric hyperbolic form as 
1 + ^ 

drip = dp (^^-^^TT - pi?j , (23) 

The initial data are arbitrary in our conjectures. In the numerical studies presented below 
we choose a Gaussian pulse 

$(0,p) = Ae-^"-"'^''/'"', 9^$(0,p) = 0, (24) 

with fixed parameters pc = 0.3, cr = 0.07, and varying amplitude A. Tests for different 
initial data give the same qualitative results which makes us feel confident that the phenomena 
described below are universal. 



4.2. The code 

We solve the hyperboloidal initial value problem ( 123124b numerically using 4th order Runge- 
Kutta integration in time and 6th order finite differencing in space. At the origin we apply the 
regularity condition ?/'(t, 0) = 0. No boundary conditions are needed at the outer boundary 
because there are no incoming characteristics due to the compactification. One-sided finite 
differencing is applied on the numerical boundaries. 

To test the code we performed a three level convergence study with 200, 400 and 800 
grid cells on the coordinate domain p e [0, 1] that has infinite physical extent. For this study, 
we used an amplitude of A = 2 which is below the critical amplitude, and a Courant factor 
of Ar/A/5 = 0.8. The convergence factors shown in Figure |2] indicate that, as expected, 
our code is 6th order convergent after a short transient phase. The dashed curve has been 
calculated using a much smaller Courant factor of 0.08 to show that the initial transient phase 
is due to numerical errors in the time integration that converge with 4th order The long 
time convergence plot on the right panel in Figure |2] shows that convergence is lost at about 
T — 1000 for the number of cells used in this study. This is due to accumulation of numerical 
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Figure 3. Left panel: The local power index calculated along r = {100,20} and 

at the origin from top to bottom. Right panel: A numerical solution that blows up at the 
origin is depicted by small squai'es and the theoretical prediction of the blowup is depicted 
by the solid line on a log-log scale. The small squai'es ai'e not distributed uniforinly in time 
because we reduce our time steps while the solution grows. The figures indicate that the code 
reproduces the known decay rates for small solutions and the blowup rate for large solutions 
very accurately. 



errors and occurs at later times in tests with higher resolution. We give numerical evidence 
for universal dynamics only in the convergent regime. Two-level convergence tests using the 
explicit conformal solution iT% give the same qualitative results. 

Further tests can be performed using known properties of solutions to the cubic wave 
equation as studied in HJ. First, it is known that generic small solutions decay as t^^ near 
timelike infinity and along null infinity It is clear from ^ that the attractor solution 
has this behavior for 5 > 0. In order to accurately measure the decay rate of the solution along 
null infinity, near timelike infinity and in the transition domain between these two asymptotic 
regimes, we calculate the local power index, Pp{t), defined by 

dlnr 

Figure [3] shows that our code reproduces the predicted decay rates very accurately. 

Second, it is known that large initial data lead to blowup in finite time with the blowup 
rate at the origin V2/{T-t) 01 H. This is in accordance with the attractor solution (|7| for 
6 < 0. The right panel in Figure [3] shows the numerical solution at the origin for large data 
on a log-log scale against the theoretical prediction of the blowup rate. The two curves match 
over many orders of magnitude indicating that our code can reliably handle the blowup. 



5. Numerical evidence for universal dynamics 



The space of solutions to the cubic wave equation can be divided into three parts: decay, 
blowup and criticality. These parts correspond respectively to 6 > 0, < 0, and = for 
the attractor solution We will follow this natural classification in our presentation of the 
numerical evidence for the universality of dynamics. 

In many cases, the numerical evidence will be presented in terms of the conformally 
rescaled solution in the coordinates presented in the previous section. In these variables the 
two-parameter family of solutions ^ takes the form (using the abbreviation f = r + a) 

" (f + (f + 1) + 1) - _ (f - 1) + 1) ■ ^^^^ 
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Figure 4. Left panel: The L2-norm of the difference between a generic decaying numerical 
solution and the fitted attractor solution is depicted by the small squares on a log-log scale. The 
solid line, shown for comparison, has the slope —4, in accordance with Conjecture 1. Right 
panel: Modulation in the relative eiTors 5 of the parameters a (solid line) and b (dashed line) 
on a log-log scale. The modulation for a seems to decay as t^^ and for b as t^^ . The relative 
error, say 5a for a, is defined as 5a = \o.{t) — a\/a where a is the value of the parameter at 
the last time step of the numerical evolution. 

5.1. Decay 

5.1.1. Convergence to the attractor According to Conjecture 1, the difference between the 
attractor solution (with optimal parameters) and a generic solution should decay as for 
small data ( fTSt . This behavior is confirmed numerically in Figure|4] 

For this plot, we first fit the numerical solution to the attractor dZSl ) to determine the 
optimal parameters (a, h). We applied two methods for the fit. Fit in time, and fit in space. 
For the fit in time we fix a grid point p =- p/, fit the solution ^{t, pf) in t and repeat the 
process for each grid point pf. This gives us a set of values (a, b)pj depending on pf. The 
average of these values over all grid points gives the optimal parameters. For this method, the 
starting time for the fit needs to be taken well behind the approximate wavefront because the 
convergence of the numerical solution to the attractor family, and therefore the variation of 
(a, b)p^ over the grid, will depend strongly on the foliation in early times. Note that our grid 
has infinite physical extent due to compactification. The variation is expected to decrease in 
time which we have checked numerically. Choosing a large value of r as the starting time of 
the time fit ensures higher accuracy in the determination of the optimal parameters. 

The second method that we applied is the fit in space in which we fix a time r = , fit 
the solution $(t/, p) in p and repeat the process on each time surface tj. This gives us a set of 
values (a, b)^^ depending on ry. The resulting time variation of the parameters is commonly 
referred to as modulation. If our conjecture is correct, the modulation should be small. 

At the end, both methods should deliver the same values up to a small error. This is used 
as a cross-check for the quality of fitting. The modulation of the parameters a and b in time is 
shown on the right panel in Figure |4] We observe that the accuracy of the fitting in b is better 
than in a. This is due to the fact that the perturbation generated by changing a decays in time 
while the perturbation generated by changing b does not (fT4l i. hence an error in determining b 
is easier to spot. 

Once the optimal parameters have been determined, we compute the difference between 
the attractor solution and the numerical solution at each time step. The difference $ — ^(a,b) 
at late times has a strict sign for all values of p on our numerical grid. Therefore the L2-norm 
of this difference over the grid as a function of time gives a strong uniform measure of how 
close the solutions are to each other. Figure|4]indicates that this difference falls off as i^^, as 
claimed in Conjecture 1. 
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Figure 5. On the left panel we plot decay rates for the flip solution with Aj = 2.4913 along 
the surfaces .y^"*", r = {100, 33, 14, 0} from top to bottom. On the right panel we plot the 
decay rates for a solution with initial amplitude Aj — 0.01 along the same surfaces. Here, the 
generic decay rates are obtained after a much longer time than that in Figure [3] 



5.7.2. Exceptional solutions As mentioned above in Remark 3, not all globally regular 
solutions converge to the attractor 0(a,h) ■ The existence of exceptional solutions with different 
asymptotic behavior can be seen as follows. Consider a one-parameter family of initial data, 
such as ( l24b . It turns out that along such a family typically there occurs a flip of sign of the 
attractor, that is, there is a (subcritical) amplitude Af, such that for A < Af the solutions 
converge to, say, 0(a.h) and for A > Af they converge to —(t>{a,b)- Performing a bisection 
we may easily fine-tune the amplitude to and generate a special solution (below referred 
to as the flip solution) that is not an element of the attractor family. One can expect that the 
flip solution will have a faster decay rate than the generic global solutions. This expectation 
is confirmed in Figure |5] which shows that the flip solution decays as t^^ at timelike infinity 
and as along null infinity. By modifying the initial amplitude slightly away from the flip 
solution, we can see that the decay rates of the generic solutions are attained at late times after 
a transient phase. Note that by construction the flip solutions correspond to codimension- 
one initial data. It is likely that there exist initial data of higher codimensions which lead to 
globally regular solutions with even faster decay rates, however the numerical construction of 
such solutions would be very difficult. 



5.2. Blowup 

5.2.1. Convergence to the attractor Merle and Zaag proved in O that the ODE solution 
y/2/{T — t) determines the universal rate of blowup for equation ([U, however the problem 
of profile of blowup remains open. Numerical simulations in spherical symmetry H showed 
that for a solution which blows up at the origin, its deviation from \/2/ {T — t) near the tip of 
the past light cone of the blowup point is very well approximated by the second eigenmode in 
the expansion (fT2b (see Figure 4 in |4| ) 

0(^,r)-^«.c,(l-^^j, (26) 

however an error on the right hand side was not quantified in IH. According to Conjecture 2 
the approximation ( |26] | may be improved to 

0(t,r)-0(„,b)(t,r)«O((T-t)2), (27) 

provided that the parameters (a, h) are optimal. 

The numerical verification of the formula (IZTT i is shown in Figure |6l The data for this 
plot were produced using a code based on standard Minkowski coordinates. 
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Figure 6. Blowup at the origin in standard coordinates on a log-log scale. The difference 
between the attractor solution and the numerical solution is depicted by small squares. The 
linear fit (solid line) gives the slope 2.02, which confirms Conjecture 2. The deviation from 
the straight line seen at times very close to the blowup time is due to the fact that the domain 
of analysis extends beyond the past light cone of the blowup point. 
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Figure 7. Simultaneous blowup along the numerical grid for fine tuned initial data. The 
solution is plotted on a log scale at various time steps close to the blowup time T. The solid 
lines depict an attractor solution with b = —1/2 at the coiTesponding times. The times are 
from bottom to top AT = T - t = {0.34, 0.1, 0.02, 0.007, 0.003, 0.0008}. 



We point out that there is no genericity condition in Conjecture 2. This might appear 
surprising in view of existence of a countable family of regular self-similar solutions of 
equation (|2|l lfT4l . It seems that these self-similar solutions do not play any role in the Cauchy 
evolution, which is probably due to the fact that they contain singularities outside the past 
light cone of the blowup point (see section 6 in 

5.2.2. Blowup surface For b < the attractor solutions blow up along a hyperboloid 

+ (28) 

This surface has the form ( l20b with the mean extrinsic curvature K = —66. Hence, for our 
choice of the hyperboloidal foliation with K = 3 the blowup surface of the attractor with 
b = —1/2 coincides with one leaf of the foliation. Therefore, if Conjecture 2 is correct, 
in the case of the blowup solution converging to the attractor with the optimal parameter 
b = —1/2 we should observe an approximately simultaneous blowup along the whole grid. 
This expectation is verified in Figure |7] 

Note that, for a given constant mean curvature foliation, the simultaneous blowup is not 
generic. To produce data for Figure |7] we used the dependence of b on the amplitude A of 
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Figure 8. Blowup at null infinity. The numerical solution (small squares) at the times 
(counting, near the origin, from top to bottom) r = {2.3, 3.2, 5, 7, 25} is compared to the 
attractor solution (solid line) with b = —0.02. We see that the solution grows at null infinity 
while it decays near the origin. 



the gaussian. For relatively small (but supercritical) values of A we have — 1/2<&<0 and 
the blowup first occurs at null infinity, while for larger amplitudes we have 6 < — 1/2 and the 
blowup first occurs at the origin. Performing bisection between these two states we fine-tuned 
to the blowup solution with b = —1/2. 

We would like to emphasize that a rather counter-intuitive phenomenon of blowup at 
null infinity is a mere coordinate effect. This only occurs if the mean extrinsic curvature 
of the blowup surface, given by —66, is smaller then the mean extrinsic curvature of the 
hyperboloidal foliation, K, used in the numerical simulation. Nevertheless, as discussed 
above, we can use this effect to our advantage to probe the shape of the blowup surface. Note 
also that one can predict the blowup at null infinity by fitting the attractor to the numerical 
solution at the origin. If this fit gives a value of 6 S (—1/2, 0), then we know that the solution 
will blow up at null infinity even though it is decaying at the origin (see Figure|8]l. 

5.3. Critical behavior 

Now we consider the behavior of solutions for initial data lying at the boundary between 
dispersion and blowup. Let us recall that this problem was addressed before in |0] for 
the focusing wave equation dtt4' ^ ^4' --4)^ — with three values of the exponent p 
(corresponding to three different criticality classes with respect to scaling of energy): p = 3 
(subcritical), p = 5 (critical), and p = 7 (supercritical). It was shown there that the 
nature of the critical solution, whose codimension-one stable manifold separates blowup from 
dispersion, depends on p: for p = 7 the critical solution is self-similar, while for p = 5 it is 
static. For p = 3 the critical solution was not determined because of numerical difficulties 
(although with hindsight it could have been inferred from Figure 1 1 in ID). Now, in view of 
Conjectures 1 and 2 which assert that the solution (t)(^a,b) is an attractor for generic dispersive 
solutions if & > and for blowup solutions if 6 < 0, it is easy to guess that the critical solution 
corresponds to 6 = 0, hence it has the form 0c = V^/{t + a). 

The critical solution is difficult to study with standard numerical methods by bisection 
because it is a globally decaying solution, but it can be studied very accurately with the 
conformal method. The reason is that the conformal scaling factors out the leading asymptotic 
behavior implying that the 1/t decay of the critical solution is factored out at the conformal 
boundary. Specifically, the rescaled critical solution $c = in the new coordinates as 

given in dZST l evaluated at null infinity becomes <l>c|j^+ = *i'(a,o)('''i 1) ~ V2, hence the 
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Figure 9. Left panel: The critical solution on the grid at times r = {2, 3, 5, 8, 13, 34, 89} 
(counting from top to bottom). We see that at r = 2 the solution is not yet described well by 
the attractor solution, but already after r = 3 the agreement is good. Right panel: Numerical 
solutions close to the critical solution evaluated at null infinity. Denoting the initial amplitude 
for the solution that corresponds to the almost constant solid line at \/2 by Ac, the deviating 
dashed curves have A = Ac it 10~®, the deviating solid ones have A = Ac it 10~*. 



deviation of $|^+ from V2 can be used as a bisection criterion. On the left panel of Figure 
|9]we plot the critical solution on the grid at various time steps and compare it to the attractor 
solution with 6 = and fitted a. We see that the solution is decaying while its value at 
^+ = {p = 1} is constant. The instability of the critical solution can be seen by evolving 
initial data that differ slightly from the critical amphtude Ac - for such data the deviation from 
V2 at grows linearly with time. This is depicted on the right panel of Figure |9] for four 
different values of marginally critical amplitude. 

6. Final remarks 

We wish to emphasize that the use of the hyperboloidal foliation (|20] | in combination with 
the conformal method was instrumental in unraveling the dynamics of global solutions 
of equation (|2]i. First and foremost, this method eliminates the need of introducing 
an artificial boundary, which is a notorious problem in computing wave propagation on 
unbounded domains. Second, the intersection of t = const hyperboloids with ^+ increases 
monotonically with r leading to the dispersive dissipation of energy along the leaves of 
the foliation which is a mechanism responsible for convergence to the attractor. Third, the 
conformal rescaling allows a very accurate computation of the critical solution by factoring 
out its leading asymptotic behavior Finally, with this approach we can probe efficiently the 
shape of the blowup surface and observe the simultaneous blowup on the whole grid. 

Our main observation that a simple family of exact solutions can act as a universal 
attractor for solutions of the nonlinear wave equation was unexpected to us. It is clear that 
this surprising phenomenon is intimately related to the conformal invariance of the cubic wave 
equation, and therefore it is more a curiosity than a stable property, in particular it is absent 
for semilinear focusing wave equations 

□(^+|0|f-V = O (29) 

with p ^ 3. However, we conjecture that the threshold behavior mediated by the 
slowly decaying global solution is structurally stable in the sense that the ODE solution of 
equation (l29T l 

1 

p-i 

(30) 



it + aV' 



where 



2(P+1) 
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is a critical solution for all exponents satisfying 1 + V2 < p < 3. Note that the decay rate 
exponent of the critical solution, 2/(1 — p), and the decay rate exponent of generic solutions, 
1 — p, merge for p ^ 1 + \/2, which is consistent with the fact that for p < 1 + \f2 all 
solutions of equation (|29T l with compactly supported initial data blow up in finite time ifTSl . 
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